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ABSTRACT 


We  present  a  new  algorithm  for  image  restoration  with  application  to  image  spectrometry,  combining 
two  radically  different  techniques:  the  singular  value  decomposition  (SVD)  and  the  method  of  projections 
onto  convex  sets  (POCS).  The  SVD  technique  is  used  to  obtain  an  initial  estimate  of  the  unknown  im¬ 
age  and  to  establish  correspondence  between  the  missing  data  and  the  spectral  description  of  the  image. 
The  iterative  method  of  convex  projections  is  then  applied  to  the  estimate,  regaining  the  missing  data 
by  enforcing  a  sequence  of  constraints  on  the  reconstructed  object.  We  report  results  of  investigations  of 
the  SVD-POCS  method  and  demonstrate  that  the  new  algorithm  leads  to  significant  improvements  in  the 
recovered  image. 


1  INTRODUCTION 

The  main  problem  of  image  reconstruction  in  infrared  image  spectrometry  is  the  missing  cone  in  the 
collected  data  cube  [1,5,6,9,16].  The  cone  evolves  around  the  zero  spatial  frequency  line,  increasing  its 
radius  for  higher  chromatic  frequencies.  The  spectrometer’s  low  spatial  frequency  ‘blindness’  results  from 
the  finite  area  of  the  focal  plane  array  and  from  the  limited  dispersion  of  the  prism,  which  prohibits 
observation  of  unbounded  objects.  The  incomplete  data  phenomenon  occurs  in  many  other  applications: 
optical  astronomy  (where  data  is  lost  due  to  atmospheric  turbulence),  electron  microscopy,  SAR  and 
medical  CAT  (where  incomplete  view  data  is  acquired).  In  either  case,  a  unique  reconstruction  of  the 
observation  object  is  impossible.  Use  of  the  Moore-Penrose  inverse  via  SVD  offers  the  minimum  norm 
least  square  (MNLS)  solution;  the  obtained  estimate  suffers,  however,  from  poor  feature  resolution  and 
contains  severe  artifacts,  resulting  from  the  high  condition  number  of  the  system  transfer  function  matrix. 

In  an  alternative  approach,  one  may  avoid  the  difficulties  inherent  in  inverting  a  singular  matrix  and 
recover  the  unknown  image  by  utilizing  known  information  about  the  image.  Typically,  an  initial  (often 
arbitrary)  guess  about  the  unknown  object  is  made  and  then  it  is  subjected  to  a  sequence  of  corrections, 
forcing  it  to  satisfy  a  number  of  desirable  characteristics,  such  as  finite  boundary,  positivity,  etc.  This 
sequence  of  corrections  is  applied  repetitively,  until  the  reconstruction  error  becomes  sufficiently  low.  At  its 
simplest,  this  idea  is  realized  in  the  form  of  the  Gerchberg-Papoulis  procedure  [14],  where  the  missing  data 
segment  is  extrapolated  in  the  Fourier  domain,  and  then  used  to  improve  the  resolution  of  the  image.  The 
basis  for  this  technique  is  that  a  spatially  limited  object  has  an  analytic  Fourier  transform,  and  therefore 
can  be  uniquely  determined  in  principle  from  any  finite  interval  of  its  spectrum.  The  Gerchberg-Papoulis 
approach  yields  an  iterative  algorithm,  which  approaches  the  solution  object  by  alternating  between  space 
domain  and  Fourier  domain.  The  known  parts  of  the  spectrum  and  the  known  boundary  of  the  image  are 
imposed  on  the  iterated  solution  as  constraints. 

This  powerful  idea  is  extended  further  by  the  method  of  projections  onto  convex  sets,  where  additional 
a  priori  knowledge  about  the  image  (non-negativity,  energy,  mean,  etc.)  is  incorporated  into  the  algorithm. 
As  a  result  an  increased  rate  of  convergence  and  an  improved  performance  is  achieved,  while  producing 
only  a  modest  increase  in  the  computational  cost  due  to  implementation  of  the  additional  constraints.  This 
advantageous  performance/complexity  trade-off  is  of  particular  importance  in  image  processing  situations, 
where  data  dimensionality  prohibits  computationally  complex  approaches.  The  main  drawback  of  the 
POCS  method  is  its  sensitivity  to  the  initial  guess.  It  is  known  that  if  the  initial  estimate  is  not  made 
carefully,  the  method  converges  very  slowly. 

We  propose  a  hybrid  algorithm  which  combines  the  SVD-based  restoration  scheme  and  the  POCS 
method  in  a  way  that  minimizes  the  shortcomings  of  both  techniques.  The  SVD  technique  is  used  to 
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obtain  an  initial  estimate  of  the  unknown  image  and  to  establish  the  correspondence  between  the  condition 
number  of  the  OTF  matrix  and  the  dimension  of  the  transform  domain  constraint.  POCS  is  then  applied 
to  the  estimate  to  reduce  blur  and  artifacts,  each  POCS  iteration  producing  a  corrected  version  of  the 
previous  estimate.  The  fusion  of  the  two  techniques  leads  to  an  efficient,  robust  and  general  image  recovery 
tool.  The  combined  SVD-POCS  approach  is  viable  whenever  the  optical  transfer  function  model  is  known 
prior  to  image  measurement  (and  therefore  SVD  computation  can  be  done  off-line).  It  yields  significant 
improvements  in  the  restored  imagery,  compared  to  the  pseudoinverse  method. 


2  PROBLEM 


All  imaging  spectrometer  reconstructs  a  three  dimensional  spatial- chromatic  object  cube  from  a  sequence  of 
two-dimensional  images.  The  reconstruction  can  be  accomplished  in  several  ways,  depending  on  whether 
multiplexing  of  the  information  is  performed  in  the  spatial  or  chromatic  domains,  or  both.  Different 
multiplexing  schemes  imply  different  trade-offs  in  terms  of  efficiency,  flexibility  and  complexity  of  the 
spectrometer;  the  numerical  formulation  of  the  reconstruction  problem  however  remains  in  either  case 
essentially  the  same. 

In  [13]  Mooney  proposed  a  direct  vision  prism  computed- tomography  image  spectrometry  technique. 
In  his  approach,  the  multiplexing  is  accomplished  by  rotating  the  prism.  As  the  prism  is  rotated,  each 
chromatic  slice  of  the  object  cube  follows  a  circular  path  with  the  radius  of  the  path  determined  by  the 
prism  dispersion.  A  sequence  of  spatial,  chromatically  overlapped  images  {^m}  is  thus  obtained,  each 
image  being  a  linear  superposition  of  all  chromatic  slices  through  the  object  cube  {/n}  spatially  convolved 
with  a  point  spread  function 

N-l 

y)=  2/)  *  *  fn{x,  y)y  (1) 

n=0 

where  x  and  y  denote  spatial  indexes.  A  Fourier  domain  equivalent  of  (1)  is 


■go(6C) 

rfo(e,c)  1 

si(6  0 

=  H(e,C) 

fi(e,C) 

.  gM-l(^>C)  - 

.fN-l(CC)- 

where  the  H(^,C)  is  an  M  x  iV  matrix 


Ho,o(?,C) 

■  Ho,n-i(CO  ■ 

Hi,o(CC) 

•  Ha,N-i(e,0 

.  HM_i,o(e,C) 

and  gm(^,C)»  Hm, n(^>C)>  fn(£)C)  represent  two-dimensional  discrete  Fourier  transforms  oi 
hm.,n{x,y)  and  fn(x,y),  respectively.  In  shorthand  (2)  takes  the  form 

g  =  Hf .  (4) 
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The  objective  of  computed- tomography  image  spectrometry  is  to  find  f,  given  the  measurements  g 
and  the  optical  transfer  function  matrix  H.  Unfortunately,  a  unique  solution  to  (4)  does  not  in  general 
exist.  First,  both  g  and  H  are  known  only  within  a  finite  resolution.  Secondly,  even  an  undistorted  g  does 
not  completely  determine  f  in  a  diffraction-limited  system.  As  a  result,  the  matrix  H  becomes  singular 
or  ill-conditioned  at  low  spatial  frequencies  and  at  high  chromatic  frequencies,  resulting  in  a  region  of 
information  ambiguity  in  the  recovered  image,  known  as  the  ‘missing  cone’  .  The  situation  becomes  even 
worse  when  the  recorded  data  g  is  contaminated  by  noise,  i.e. 

g  =  Hf  +  n.  (5) 

One  possible  approach  which  allows  exercising  some  degree  of  control  over  instabilities  arising  from 
inverting  an  iU- conditioned  matrix,  is  application  of  the  singular  value  decomposition  (SVD)  to  H 

H  =  UwV+,  •  (6) 

where  U  is  an  M  x  N  column-orthonormal  matrix,  w  is  an  N  x  N  diagonal  matrix  and  V  is  an  iV  x  iV 
orthonormal  matrix.  Since  the  singularity  of  H  manifests  itself  in  zero  singular  values  of  w,  we  can  easily 
compute  a  pseudoinverse  of  H  as 

=  Vw*U+, 

where  an  element  w*  of  the  diagonal  matrix  w*  is  defined  as 

1/w,  ts  /  0, 

0,  u;  =  0. 


(7) 

(8) 


Applying  (7)  to  the  noisy  situation  (5)  allows  computing  an  estimate  of  f 

(9) 
(10) 


f  =  Vw*Utg 

=  Vw*(wVtf +  Utn). 


The  drawback  of  the  pseudoinverse  approach  is  that  if  one  or  more  entries  in  the  diagonal  matrix  w 
is  zero,  then  there  is  insufficient  information  to  determine  how  the  measured  data  should  be  distributed 
among  chromatic  bands  and  as  a  result  an  intra-band  crosstalk  occurs.  Additionally,  for  small  values  in  w, 
the  noise  in  (10)  is  going  to  be  greatly  amplified.  One  way  to  reduce  the  noise  amplification,  is  to  modify 
(8) 


1/w,  W  >  €thr 
0,  W  <I  f-thri 


(11) 


where  ethr  is  a  threshold,  so  that  the  trade-off  between  the  crosstalk  and  noise  amplification  can  be 
balanced.  Alternatively,  one  can  apply  a  form  of  Wiener  filtering  to  (11),  replacing  w*  with  w/(w'^ + 

In  either  case  the  balancing  between  the  crosstalk  and  noise  amplification  is  obstructed  by  the  fact  that 
the  threshold  can  only  be  roughly  estimated,  as  the  image/noise  statistics  are  rarely  known  with  sufficient 
precision. 

Two  regularization  approaches  to  the  iU-posed  (5)  are  possible:  one  is  based  on  statistical  models  of 
the  data  (i.e.  maximum  likelihood,  maximum  entropy,  etc.  [10,20]),  and  another  is  deterministic,  utilizing 
available  a  priori  information  about  the  unknown  image.  In  the  next  section,  we  will  introduce  an  iterative 
algorithm  of  projections  onto  convex  sets,  belonging  to  the  second  class  of  regularization  techniques. 
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3  POCS 


The  method  of  projections  onto  convex  sets  is  a  recursive  algorithm  for  finding  a  point  in  the  intersection 
of  a  given  sequence  of  closed  convex  sets.  The  method  was  developed  by  Bregman  [3]  and  Gubin  [8]  and 
adapted  to  signal  processing  by  Youla  [18,19],  POCS  was  applied  by  Levi  and  Stark  to  image  restoration 
from  phase  [11],  to  image  restoration  from  magnitude  [12],  by  Sezan  and  Stark  to  limited-angle  data 
restoration  [15]  and  by  Trussell  to  X-ray  fluorescence  signals  [17], 

In  the  simplest  case  (as  will  be  seen  later)  POCS  reduces  to  the  well-known  Gerchberg-Papoulis  method. 
We  will  provide  a  summary  of  the  POCS  method  and  list  the  relevant  convex  sets  without  proofs.  For 
details  the  reader  is  directed  to  [19], 

The  objective  of  POCS  is  to  find  an  unknown  function  f  restricted  to  lie  in  the  intersection  of  a  given 
sequence  of  R  closed  convex  sets 

Co  =  n  Cr.  (12) 

r=l 

A  subset  C  of  Ti,  where  H  is  a  Hilbert  space,  is  convex,  if  for  any  two  of  its  elements  xi  and  X2  it  contains 
the  element  x  ~  fixi  +  (1  —  iJ^)x2y  where  0  ^  £  1,  A  subset  C  of  77  is  closed,  if  the  limit  point  of  anj 

sequence  of  points  in  C  is  contained  in  C.  Given  projection  operators  Pr  associated  with  closed  convex 
spaces  Cr  a  sequence  of  functions  {fk}  is  generated  by  the  recursive  relation 

R 

fk+l  -  PRPR-l^-^Plfk  =  (11  Pr)fk.  (13) 

r=l 

such  that 

,lim  fk  =  /', 

where  is  an  element  in  Co  that  is  closest  to  /.  If  Co  contains  a  single  point,  the  solution  is  unique 
contains  more  than  a  single  point,  POCS  yields  a  feasible  solution. 

The  relation  (13)  can  be  stated  more  generally 

fk+l=if[Tr)fk,  (15) 

r=l 

where  Tr  =  I  +  XAPr  - 1),  0  <  <  2  and  /  is  the  identity  operator.  The  A,.’s  are  relaxation  parameters, 

and  can  be  used  to  accelerate  the  rate  of  convergence  of  the  algorithm.  In  general,  identification  of  opti¬ 
mal  relaxation  parameters,  i.e.  parameters  providing  the  maximum  acceleration,  is  difficult,  and  usually 
suboptimal  heuristic  choices  are  made.  In  a  simple  situation,  however,  where  only  two  closed  convex  sets 
are  used,  and  one  of  the  projections  is  linear,  a  near  optimal  per-cycle  relaxation  parameter  for  the  second 
projection  can  be  calculated  [12],  that  does  not  require  knowledge  of  the  original  data  /. 


(14) 
If  it 
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The  following  projection  operators  are  typically  used: 

1.  Ci:  The  set  of  all  functions  in  Ti  that  vanish  outside  a  region  B  C  R? 


Pif  =  XBi>J)f(iJ)  = 


/(*,/),  (*,;■)€  i?, 

0, 


(16) 


where  B^  is  a  complement  of  B  and 


Xsii  i)  -  I  ^ 


(17) 


2.  C2:  The  set  of  all  functions  in  H  whose  Fourier  transform  vanishes  outside  a  frequency  region 
AcB? 


P2f 


F{u,v),  iu,v)£A, 


0,  (ii,  u)  e 

In  the  particular  case,  when  only  Pi  and  P2  are  used,  i.e. 

fk+i  =  -P2-P1/*, 


(18) 


(19) 


POCS  degenerates  to  the  well-known  Gerchberg-Papoulis  algorithm.  Equation  (13)  can  then  be  rewritten 
as 

A+i  =  F-^Fo  +  (1  -  XA)F{xBfk}},  (20) 

where  Fq  is  the  known  part  of  the  spectrum  of  /,  and  /o  =  is  the  initial  estimate  of  /. 


3.  C3:  The  set  of  all  non-negative  functions  in 


P3f  = 


f,  />o. 


0,  /<0. 

4.  C4:  The  set  of  aU  functions  in  77,  whose  amplitudes  lie  in  a  closed  interval  [a,  6],  0  <  a  <  6 


(21) 


^4/=  ^ 


fa,  f  <  a, 
f,  a<  f  <b, 
b,  f>b. 


(22) 


Notice  that  C3  is  a  subset  of  C4. 

The  following  projections  were  additionally  considered  in  this  work: 


5.  Cjr:  The  set  of  all  non-negative  area-scaled  functions  in  H 

Ppf  =  I  /  >  0) 

\  0,  /<0, 

where  S  =  f  f  f(x,  y)dxdy  is  the  area  of  /  auid  S'^  is  the  area  of  P3/. 


(23) 
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We  have  selected  the  constraints  Pi,  P4,  Pa  and  Pr  or  Pt>  to  be  included  in  the  image  recovery 
algorithm.  The  finite  support  constraint  Pi  is  realized  by  projecting  an  X  x  Y  image  through  an  X  x 
Y'  (X'  <  X,  Y'  <  Y)  field  stop.  At  each  iteration,  pixels  residing  outside  the  field  stop  boundary  are 
reset  to  a  fixed  precomputed  level.  The  Pi  constraint  can  be  combined  with  the  Pa  constraint,  which 
corrects  the  zero  spatial  frequency  of  the  iterated  image.  The  actual  value  of  the  dc  level  of  an  image  is 
obtained  during  camera  calibration.  Additionally,  the  P4  constraint  can  be  included  in  the  group  of  space 
domain  constraints,  if  estimates  of  image  magnitude  limits  are  available.  ^ 

In  the  transform  domain  we  have  a  choice  between  Py  and  Pt'.  Pt'  is  similar  to  P2,  except  that  the 
spectral  correction  takes  place  in  the  space  spanned  by  columns  of  the  V  matrix.  The  values  of  the  product 
Vifo  are  updated  for  all  triplets  (^,  C,  n)  contained  in  the  ’missing  cone’  region  of  the  transform  space, 
defined  by  the  condition  w  <  ethr-  Pt  is  a  version  of  Py/,  which  allows  for  smooth  transition  between 
bands  (we  found  that  the  Py  constraint  performs  significantly  better).  ^  •  i 

The  image  reconstruction  algorithm  consists  of  two  stages  (Fig.l.):  first,  computation  of  the  initial 
estimate;  second,  obtaining  a  refinement  of  the  solution  by  a  series  of  POCS  itera,tions,  in  which  the 
recovered  image  is  forced  to  satisfy  predefined  characteristics.  In  the  SVD  stage,  the  image  estimate  f©  is 
obtained  by  taking  a  product  of  the  adjoint  of  matrix  H  and  the  coUected  data  vector  g,  with  rows  of  U^g 
corresponding  to  small  singular  values  being  attenuated  by  the  factor  +  e^).  Matrix  V  and  vector 

V+fo  are  stored  for  later  use  in  the  POCS  stage.  An  inverse  2-D  DFT  is  applied  to  fo,  yielding  a  starting 

vector  for  the  iteration.  n-um 

The  POCS  iteration  includes  a  sequence  of  space  domain  projections  followed  by  a  forward  2D  iih  i ,  a 
transform  domain  projection,  and  an  inverse  2D  DFT.  The  sequence  of  space  domain  projections  include 
Pi,  P4  and  Pa-  The  transform  domain  projection  is  either  Pt'  or  its  ’soft’  version  Py.  This  stage  requires 
two  operations  of  matrix-vector  multiplications:  one  for  obtaining  the  product  and  one  for  computing 
the  estimate  fk  =  V(V^fl{;). 

The  multiplicative  complexity  of  the  iterative  part  of  the  algorithm  is 

2kXYN{N  -b  1  +  log-iiXY))  (27) 

where  k  is  the  number  of  iterations,  N  is  the  number  of  colors,  and  X  X  Y  is  the  image  size. 
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V^fo 


Figure  1:  Flowchart  of  the  SVD-POCS  algorithm. 


4  RESULTS 

We  have  tested  two  sets  of  synthesized  images:  the  ‘checkerboard’  and  the  ‘point  source’.  The  size  of  both 
synthetic  image  arrays  was  NxXxY  =  8x  240  x  240,  with  the  field  stop  being  64  x  64.  The  ‘checkerboard’ 
image  contained  a  pattern  of  three  sets  of  pixels  of  intensities:  10,  300  and  1000.  The  ‘point  source’  image 
had  a  uniform  background  plane  with  the  set  of  color  intensities  {650,600,550,500,450,400,350,300},  and 
a  two-pixel  target  with  the  set  of  color  intensities  {650,600,550,500,450,400,1000,300}.  In  both  cases  10 
units  of  white  gaussian  noise  was  added  to  the  test  image. 

We  have  tested  two  implementations  of  the  POCS  algorithm 

T1  :  fk+i  =  TrTifk,  =  1, 1.9 


and 


TAF  :  fk+i  =  TxTATpfk,  Xt,a,f  =  1, 1-9 

The  basis  for  comparing  results  for  the  synthesized  images  was  the  rms  error,  defined  as 

.  HA -/II 
li/ll  ■ 

In  the  case  of  the  ‘point  source’  image,  the  T1  algorithm  {Xl  =  1)  yielded  63%  reduction  in  the  rms 
error  after  5  iterations  and  84%  after  10  iterations  (Table  1.).  Setting  the  relaxation  parameter  Ax,i  to 
1.9  increased  the  rate  of  convergence  of  the  algorithm  approximately  by  a  factor  of  two.  The  accelerated 
algorithm  started  to  diverge  at  a  moderate  rate  after  about  ten  iterations. 

The  TAF  algorithm  produced  spectacularly  good  results  yielding  89%  improvement  by  the  second 
iteration,  93%  by  the  fifth  (notice  narrowing  of  the  background  intensity  spread,  reduction  of  blocking 
effects  around  the  boundary  and  a  large  adjustment  of  the  mean  in  all  frames  in  Fig.  2.),  gracefully 
diverging  thereafter.  Increasing  the  relaxation  parameter  Xtaf  to  1.9  actually  decreased  the  rate  of 
convergence,  although  the  results  were  stiU  better  than  the  results  of  the  T1  algorithm. 

In  case  of  the  highly  textured  ‘checkerboard’  image  the  improvement  yielded  by  both  iterative  algo¬ 
rithms  was  less  radical:  45%  after  5  iterations  and  54%  after  10  iterations  for  T1  algorithm  with  =  1.0 
(Table  2.),  resulting  in  better  separation  of  the  three  groups  of  pixels  (Fig.  3.).  Similarily,  as  in  the  case 
of  the  ‘point  source’  image,  setting  Ay,!  =  1.9  provided  an  approximately  twofold  increase  in  the  rate 
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Table  1:  Summary  of  results  for  image  restoration  algorithms  T1  and  TAP  (the  ‘point  source’) 


Iteration 

number 

T1 

TAF 

Ar,i  =  1 

~  1*9 

Ar,A,F  -  1 

^T,A,F  =1.9 

0 

0.1885 

0.1885 

0.1885 

0.1885 

1 

0.1495 

0.1142 

0.0392 

0.0542 

2 

0.1225 

0.0655 

0.0214 

0.0227 

3 

0.1009 

0.0407 

0.0158 

0.0169 

4 

0.0834 

0.0285 

0.0143 

0.0163 

5 

0.0693 

0.0232 

0.0140 

0.0163 

10 

0.0308 

0.0193 

0.0152 

0.0184 

15 

0.0205 

0.0209 

0.0168 

0.0208 

20 

0.0194 

0.0227 

0.0182 

0.0228 

25 

0.0201 

0.0244 

0.0195 

0.0246 

of  convergence.  Somewhat  surprisingly  the  results  of  TAF  and  T1  algorithms  were  nearly  identical;  the 
additional  constraint  did  not  help  nearly  as  much  as  for  the  point  source  image.  f  4.^  ti 

In  short,  the  experiments  confirmed  that  increasing  the  relaxation  parameters  m  the  case  of  the  i  i 
algorithm  improve  the  rate  of  convergence  significantly.  Computation  of  the  optimum-per- cycle  parameter 
as  given  in  [12]  yielded  marginal  improvement  in  a  noisy  situation,  although  earlier  experiments  have  shown 
excellent  results  in  a  noise-free  situation.  Further  research  in  this  area  is  needed. 

Results  of  the  TAF  algorithm  have  shown  that  certain  constraints  can  yield  dramatic  improvement  lor 
one  type  of  imagery,  while  having  virtually  no  effect  on  others  (in  our  case  the  constrmnt  on  zero  spati 
frequency  had  a  convincing  effect  on  an  image  with  low  spatial  frequency  characteristics).  Even  though  the 
‘typical’  image  might  be  complex  in  character,  design  of  an  adaptive  algorithm  with  a  set  of  adjustable, 
scenery-dependent  constraints  might  be  a  reasonable  course  of  action.  In  absence  of  information  about 
the  image  spectrum  envelope,  it  is  prudent  to  use  a  full  set  of  constraints  and  properly  adjusted  relaxation 

parameters. 


5  SUMMARY 

In  this  work,  we  proposed  a  hybrid  image  restoration  algorithm  which  combines  the  direct  inverse  approach 
with  an  iterative  technique  that  uses  a  priori  information  about  the  image.  We  have  tested  two  iterative 
methods:  the  acceleration  of  the  Gerchberg-Papoulis  algorithm  with  variable  relaxation  parameters  and 
the  POCS  algorithm  in  which  an  additional  a  priori  information  about  low  spatial  frequencies  01  the 
data  vector  was  utilized.  Early  simulations  have  shown  that  application  of  the  hybrid  algorithm  can  ea 
to  significant  improvements  in  data  restoration  results.  It  was  demonstrated  that  judicious  selection  of 
relaxation  parameters  can  impact  the  rate  of  convergence  of  the  two  iterative  techniques,  even  in  the 
presence  of  moderate  levels  of  noise.  The  theory  of  convex  projections  provided  a  convenient  framework  by 
which  additional  a  priori  information  about  the  data  vector  could  be  used  to  obtain  a  high  quality  image 
in  a  limited  number  of  iterations. 
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Table  2:  Summary  of  results  for  image  restoration  algorithms  T1  and  TAF  (the  ‘checkerboard’). 


Iteration 

number 

T1 

TAF 

Ar,i  =  1 

—  1.9 

^T,A,F  =  1 

^T,A,F  =1-9 

0 

0.2245 

0.2245 

0.2245 

0.2245 

1 

0.1796 

0.1538 

0.1794 

0.1539 

2 

0.1562 

0.1295 

0.1560 

0.1297 

3 

0.1416 

0.1147 

0.1415 

0.1147 

4 

0.1315 

0.1059 

0.1313 

0.1055 

5 

0.1240 

0.1000 

0.1238 

0.0995 

10 

0.1037 

0.0896 

0.1036 

0.0895 

15 

0.0961 

0.0916 

0.0960 

0.0921 

20 

0.0938 

0.0966 

0.0938 

0.0975 

25 

0.0943 

0.1019 

0.0943 

0.1029 

While  the  SVD-POCS  approach  shows  promise,  further  work  in  tuning  performance  of  the  algorithm 
is  obviously  needed.  One  parameter,  that  was  observed  to  have  a  large  elfect  on  the  rate  of  convergence 
of  the  iterations,  was  the  singular  value  threshold.  Future  research  should  include  design  of  more  sophis¬ 
ticated  criteria  for  threshold  selection,  including  computation  of  a  spatially-weighted  optimum-per-cycle 
threshold.  Another  possibility  is  implementation  of  a  multiresolution  scheme,  in  which  instead  of  using  an 
optimum  threshold,  a  weighted  sum  of  images  computed  at  different  threshold  levels  is  taken  [7].  Finally, 
computational  savings,  as  well  as  error  reduction,  can  be  potentially  achieved  by  bringing  the  iterations  to 
a  joint  time-frequency  space,  thereby  gaining  control  over  the  uncertainty  principle  [4].  These  issues  are 
currently  under  investigation. 
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Figure  2:  Illustration  of  the  POCS  algorithms  performance  for  the  ‘point  source’  image  (the  seventh  frame 
of  data  and  its  histogram),  (a)  result  of  the  SVD  algorithm,  (b)  result  of  the  T1  algorithm  for  A  =  1.9  (5th 
iteration),  (c)  result  of  the  TAF  algorithm  for  A  =  1  (5th  iteration),  (d)  comparison  of  the  chromatic  series 
of  the  central  pixel  for  the  SVD  algorithm  (dotted  line),  T1  algorithm  (dashed  line)  and  TAF  algorithm 
(dashed-dotted  line).  Solid  line  marks  intensity  levels  of  the  original  noise-free  image. 
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Figure  3:  Illustration  of  the  POCS  algorithms  performance  for  the  'checkerboard’  image  (the  first  frame 
of  data  and  its  histogram),  (a)  result  of  the  SVD  algorithm,  (b)  result  of  the  T1  algorithm  for  A  =  1.9 
(10th  iteration),  (c)  result  of  the  TAF  algorithm  for  A  =  1.9  (10th  iteration),  (d)  comparison  of  the 
chromatic  series  of  the  central  pixel  for  the  SVD  algorithm  (dotted  Une),  T1  algorithm  (dashed  line)  and 
TAF  algorithm  (dashed-dotted  line).  Solid  line  marks  intensity  levels  of  the  original  noise-free  image. 


